PREPRINT 



Incorporating Forcing Terms in Cascaded Lattice-Boltzmann 
Approach by Method of Central Moments 

Kannan N. PremnathQ 

Department of Chemical Engineering, 
University of California, Santa Barbara, 
Santa Barbara, CA 93106 
and 

MetaHeuristics LLC, 3944 State Street, 
Suite 350, Santa Barbara, CA 93105 

Sanjoy BanerjeC 
Department of Chemical Engineering 
Department of Mechanical Engineering 
Bren School of Environmental Science and Management 
University of California, Santa Barbara, 
Santa Barbara, CA 93106 



(Dated: February 29, 2012) 



1 



Abstract 

Cascaded lattice-Boltzmann method (Cascaded-LBM) employs a new class of collision operators 
aiming to stabilize computations and remove certain modeling artifacts for simulation of fluid flow 
on lattice grids with sizes arbitrarily larger than the smallest physical dissipation length scale (Geier 
et al, Phys. Rev. E 63, 066705 (2006)). It achieves this and distinguishes from other collision 
operators, such as in the standard single or multiple relaxation time approaches, by performing 
relaxation process due to collisions in terms of moments shifted by the local hydrodynamic fluid 
velocity, i.e. central moments, in an ascending order-by-order at different relaxation rates. In this 
paper, we propose and derive source terms in the Cascaded-LBM to represent the effect of external 
or internal forces on the dynamics of fluid motion. This is essentially achieved by matching the 
continuous form of the central moments of the source or forcing terms with its discrete version. 
Different forms of continuous central moments of sources, including one that is obtained from a 
local Maxwellian, are considered in this regard. As a result, the forcing terms obtained in this new 
formulation are Galilean invariant by construction. To alleviate lattice artifacts due to forcing terms 
in the emergent macroscopic fluid equations, they are proposed as temporally semi-implicit and 
second-order, and the implicitness is subsequently effectively removed by means of a transformation 
to facilitate computation. It is shown that the impressed force field influences the cascaded collision 
process in the evolution of the transformed distribution function. The method of central moments 
along with the associated orthogonal properties of the moment basis completely determines the 
analytical expressions for the source terms as a function of the force and macroscopic velocity 
fields. In contrast to the existing forcing schemes, it is found that they involve higher order terms 
in velocity space. It is shown that the proposed approach implies "generalization" of both local 
equilibrium and source terms in the usual lattice frame of reference, which depend on the ratio of 
the relaxation times of moments of different orders. An analysis by means of the Chapman-Enskog 
multiscale expansion shows that the Cascaded-LBM with forcing terms is consistent with the 
Navier-Stokes equations. Computational experiments with canonical problems involving different 
types of forces demonstrate its accuracy. 

PACS numbers: 47.11. Qr,05.20.Dd,47.27.-i 
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I. INTRODUCTION 



Lattice-Boltzmann method (LBM), based on minimal discrete kinetic models, has at- 
tracted considerable attention as an alternative computational approach for fluid mechanics 
problems Q-Q. While its origins can be traced to lattice gas automata [5[ as a means to re- 
move its statistical noise [6], over the years, the LBM has undergone major series of advances 
to improve its underlying models for better physical fidelity and computational efficiency 
Moreover, its connection to the continuous Boltzmann equation as a dramatically simplified 
version [7|, |8[ established it as an efficient approach in computational kinetic theory and led 
to the development of asymptotic tools 9j providing a rigorous framework for numerical 
consistency analysis. The LBM is based on performing stream-and-collide steps to compute 
the evolution of the distribution of particle populations, such that its averaged behavior 
recovers the dynamics of fluid motion. The streaming step is a free-flight process along 
discrete characteristic particle directions designed from symmetry considerations, while the 
collision step is generally represented as a relaxation process of the distribution function to 
its attractors, i.e. local equilibrium states. Considerable effort has been made in developing 
models to account for various aspects of the collision process, as it has paramount influence 
on the physical fidelity and numerical stability of the LBM. 

One of the simplest and among the most common is the single-relaxation-time (SRT) 
model proposed by Chen et al. 10] and Qian et al. which is based on the BGK ap- 
proximation On the other hand, d'Humieres (1992) proposed a moment method, 
in which various moments that are integral properties of distribution functions weighted by 
the Cartesian components of discrete particle velocities of various orders are relaxed to their 
equilibrium states at different rates during collision step, leading to the multiple-relaxation- 
time (MRT) model. It is an important extension of the relaxation LBM proposed earlier 



by Higuera et al 
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15] . While it is a much simplified version of the latter, the major 



innovation lies in representing the collision process in moment space |16[ rather than the 
usual particle velocity space. By carefully separating the relaxation times of hydrodynamic 
and non-hydrodynamic moments, it has been shown that the MRT-LBM significantly im- 
proves the numerical stabi.it y HQ and better physica. representation in certain problems 
such as kinetic layers near boundaries 19], when compared with the SRT-LBM. Such MRT 



models have recently been shown to reproduce challenging fluid mechanics problems such 
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as complex turbulent flows with good quantitative accuracy 20|, l2l). An important and 



natural simplification of the MRT model is the two-relaxation-time (TRT) model, in which 



the moments of even and odd orders are relaxed at different rates 
From a different perspective, Karlin and co-workers 



221. 



23H27| have developed the so-called 



entropic LBM in which the collision process is modeled by assuming that distribution func- 
tions are drawn towards their attractors, which are obtained by the minimization of a 
Lyapanov-type functional, i.e. the so-called H-theorem is enforced locally, while modulating 
the relaxation process with a single relaxation time to maintain numerical stability. It may 
be noted that in contrast to the SRT or MRT collision operators, which employ equilibria 
that are polynomials in hydrodynamic fields, the entropic collision operator, in general, re- 
quires use of non-polynomial or transcendental functions of hydrodynamic fields. Recently, 
using this framework, a novel entropy-based MRT model was derived [28] and a Galilean 
invariance restoration approach was developed 29|. In addition, there has been consider- 



able progress in the development of systematic procedures for high-order lattice-Boltzmann 



models 
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Recently, Geier et al. 32[ introduced another novel class of collision operator leading to 
the so-called Cascaded-LBM. Collision operators, such as the standard SRT or MRT models, 
are generally constructed to recover the Navier-Stokes equations (NSE), with errors that are 
quadratic in fluid velocity. Such models, which are Galilean invariant up to a lower degree, 
i.e., the square of Mach number, are prone to numerical instability, which can be alleviated 
to a degree with the use of the latter model. Recognizing that insufficient level of Galilean 
invariance is one of the main sources of numerical instability, Geier proposed to perform 
collision process in a frame of reference shifted by the macroscopic fluid velocity. Unlike other 
collision operators which perform relaxation in a special rest or lattice frame of reference, 
Cascaded-LBM chooses an intrinsic frame of reference obtained from the properties of the 
system itself. The local hydrodynamic velocity, which is the first moment of the distribution 
functions, is the center of mass in the space of moments. A coordinate system moving locally 
with this velocity at each node is a natural framework to describe the physics of collisions 
in the space of moments. This could enable achieving a higher degree of Galilean invariance 
than possible with the prior approaches. It may be noted that the moments displaced by 
the local hydrodynamic velocity are termed as the central moments and are computed in a 
moving frame of reference. On the other hand, the moments with no such shift are called 
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the raw moments, which are computed in a rest frame of reference. 

Based on this insight, the collision operator is constructed in such a way that each central 
moment can be relaxed independently with generally different relaxation rates. However, it 
is computationally easier to perform operations in terms of raw moments. Both forms of 
moments can be related to one another in terms of the binomial theorem, and hence the 
latter plays an important role in the construction of an operational collision step. As a result 
of this theorem, central moment of a given order are algebraic combinations of raw moments 
of different orders, with their highest order being equal to that of the central moment. In 
effect, the evolution of lower order raw moments influences higher order central moments and 
not vice versa. Thus, due to this specific directionality of coupling between different central 
and raw moments, starting from the lowest central moment, we can relax successively higher 
order central moments towards their equilibrium, which are implicitly carried out in terms of 
raw moments. Such structured sequential computation of relaxation in an ascending order 
of moments leads to a novel cascaded collision operator, in which the post-collision moments 
depend not only on the conserved moments, but also on the non-conserved moments and on 
each other. 

Moreover, it was found that relaxing different central moments differently, certain arti- 
facts such as aliasing that cause numerical instability for computation on coarse grids, whose 
sizes can be arbitrarily larger than the smallest physical or viscous dissipation length scale 
can be avoided. In particular, this is achieved by setting the third-order central moments 
to its equilibrium value, while allowing only the second-order moments to undergo over- 



relaxtion 
condition 



33 



The limit of stability is now dictated only by the Courant-Friedrichs-Lewy 
34] typical of explicit schemes and not by effects arising due to the discreteness of 
the particle velocity set. Prevention of such ultra-violet catastrophe in under-resolved com- 
putations could enable application of the LBM for high Reynolds number flows or for fluid 
with low viscosities. Further insight into the nature of the gain in numerical stability with 
Cascaded-LBM is achieved with the recognition that unlike other collision operators which 
appear to introduce de-stabilizing negative hyper-viscosity effects that are of second-order 
in Mach number due to insufficient Galilean invariance, the former seems to have stabilizing 
positive and smaller hyper- viscosity effects that are of fourth-order in Mach number 35 ]. 



Recently, Asinari 



36| showed that cascaded relaxation using multiple relaxation times is 



equivalent to performing relaxation to a "generalized" local equilibrium in the rest frame 



6 



of reference. Such generalized local equilibrium is dependent on non-conserved moments as 
well as the ratio of various relaxation times. 

Clearly, several situations exist in which the dynamics of fluid motion is driven or affected 
by the presence of external or self-consistent internal forces. Examples include gravity, 
magnetohydrodynamic forces, self-consistent internal forces in multi-phase or multi-fluid 
systems. Moreover, subgrid scale (SGS) models for turbulence simulation can be explicitly 
introduced as body forces in kinetic approaches [21], Q • Thus, it is important to develop 
a consistent approach to introduce the effect of forces that act on the fluid flow in the 



Cascaded-LBM. 
for example, in 



he method for introducing force terms in other LBM approaches are given, 
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41], in which notably Guo et al. j4l| developed a consistent approach 



which avoided spurious effects in the macroscopic equations resulting from the finiteness of 
the lattice set. 

The approach proposed in this paper consists as follows. It consists of deriving forcing 
terms which can be obtained by matching their discrete central moments to their corre- 
sponding continuous version. In this regard, we consider two different sets of ansatz for 
the continuous source central moments - one based on a continuous local Maxwellian and 
another one which makes specific assumptions regarding the effect of forces for higher order 
moments. An important feature of our approach is that by construction the source terms are 
Galilean invariant, which would be a very desirable aspect from both physical and computa- 
tional points of view. To facilitate computation, the central source moments are related to 
corresponding raw moments, which are, in turn, expressed in velocity space. Furthermore, 
to improve temporal accuracy, the source terms are treated semi-implicitly. The implicit- 
ness, then, is effectively removed by applying a transformation to the distribution function. 
A detailed a priori derivation of this central moment method is given so that it provides 
a mathematical framework which could also be useful for extension to other problems. We 
then establish the consistency of our approach to macroscopic fluid dynamical equations by 
performing a Chapman- Enskog multiscale moment expansion. It will be shown that when 
Cascaded-LBM with forcing terms is reinterpreted in terms of the rest frame of reference (as 
usual with other LBM), it implies considering a generalized local equilibrium and sources, 
which also depend on the ratio of the relaxation times of various moments, for their higher 
order moments. Numerical experiments will also be performed to confirm the accuracy of 
our approach for flows with different types of forces, where analytical solutions are available. 
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This paper is structured as follows. Section |TT] briefly discusses the choice of moment basis 
employed in this paper. In Sec. IHH continuous forms of central moments for equilibrium and 
sources (for a specific ansatz) are introduced. The Cascaded-LBE with forcing terms are 
presented in Sec. IIVI In Sec. |V] we discuss the details of an analysis and the construction of 
the Cascaded-LBM and the analytical expressions for source terms. Section IVT1 provides the 
details of how the computational procedure is modified with the use of a different form of 
the central source moments. The computational procedure for Cascaded-LBM with forcing 
is provided in Sec. I VIII Results of the computational procedure for some canonical prob- 
lems are presented in Sec. IVIII1 Summary and conclusions of this work are described in 
Sec. IIXI Consistency analysis of the central moment method with forcing terms by means of 
a Chapman- Enskog multiscale moment expansion is presented in Appendix [A] Appendix [Bl 
shows that Cascaded-LBM with forcing terms is equivalent to considering a generalized local 
equilibrium and sources in the rest frame of reference. Finally, Appendix ICl investigates the 
possibility of introducing time-implicitness in the cascaded collision operator. 



II. CHOICE OF BASIS VECTORS FOR MOMENTS 



For concreteness, without losing generality, we consider, the two-dimensional, nine ve- 
locity (D2Q9) model, which is shown in Fig. [TJ The particle velocity ~£ a may be written 

as 



' (0,0) a = 



et = < 



(±1,0),(0,±1) a = l,...,4 (1) 
(±1,±1) a = 5, •••,8 

Here and henceforth, we employ Greek and Latin subscripts for particle velocity directions 
and Cartesian coordinate directions, respectively. Moments in the LBM are discrete integral 
properties of the distribution function f a , i.e. Yla=o e ax e ayfa, where m and n are integers. 
Since the theory of the moment method draws heavily upon the associated orthogonality 
properties, for convenience, we employ the Dirac's bra-ket notation in this paper. That is, 
we denote the "bra" operator (0| to represent a row vector of any state variable <fi along each 
of the particle directions, i.e. (0o, <t>i, 4>2, ■ ■ ■ , <f>s), and the "ket" operator \<f>) represents a 
column vector, i.e. (</> , <fii, 02> • • • > 4>sY , where the superscript "f" is the transpose operator. 
In this notation, (<p\(p) represents the inner-product, i.e. Xla=o ^a^a- To obtain a moment 
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FIG. 1: Two-dimensional, nine-velocity (D2Q9) Lattice. 

space of the distribution functions, we start with a set of the following nine non- orthogonal 
basis vectors obtained from the combinations of the monomials e^e™ in an ascending order. 



ip) = ir^ Q i°) = (1,1,1,1,1,1,1,1,1^, (2) 

\e ax ) = (0,1,0,-1,0,1,-1,-1,1)*, (3) 

\e ay ) = (0,0,1,0,-1,1,1,-1,-1)*, (4) 

kL + O = (0,1,1,1,1,2,2,2,2)*, (5) 

IC-O = (0,1,-1,1,-1,0,0,0,0)*, (6) 

\e ax e ay ) = (0,0,0,0,0,1,-1,1,-1)*, (7) 

K x e ay ) = (0,0,0,0,0,1,1,-1,-1)*, (8) 

\e ax e 2 ay ) = (0,0,0,0,0,1,-1,-1,1)*, (9) 

\e 2 ax e 2 ay ) = (0,0,0,0,0,1,1,1,1)*. (10) 



To facilitate analysis, the above set of basis vectors is transformed into an equivalent 
orthogonal set of basis vectors by means of the standard Gram-Schmidt procedure in the 
increasing order of the monomials of the products of the Cartesian components of the particle 
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velocities: 



\Ko) = \p) , 

\Kx) = \e ax ) , 

\K 2 ) = \e ay ) , 

\Ks) = 3|eL + e a ,)-4|p), 

\Ka\ = \p 2 -p 2 ) 

\ K &) = -3 \e 2 ax e ay ) + 2 \e ay ) , 

1-^7) = — 3 le^e^) + 2 |e ax ) , 

l^s) = 9 \e 2 ax e 2 ay ) - 6 \e 2 QX + e 2 ay ) + 4 \p) 



(11) 
(12) 

(13) 
(14) 
(15) 
(16) 
(17) 
(18) 
(19) 



This is very similar to that used by Geier et al. [321 ]. except for the negative sign used in 
\K§) by the latter. The purpose of using a slightly different orthogonal basis than that 
considered in [32| is simply to illustrate how it changes the details of the cascaded collision 
operator. It is obvious that we can define different sets of orthogonal basis vectors that differ 
from one another by a constant factor or a sign. Furthermore, it is noteworthy to compare 
the ordering of basis vectors used for the central moment method with that considered by 
Lallemand and Luo 17| : Here, the ordering is based on the ascending powers of moments 



(i.e. zeroth order moment, first order moments, second order moments,. . .) while [l7j order 
their basis vectors based on the character of moments, i.e. increasing powers of their tensorial 
orders (i.e. scalars, vectors, tensors of different ranks,. . .). 

The orthogonal set of basis vectors can be written in terms of the following matrix 



K = [\K ) , \K X ) , \K 2 ) , \K 3 ) , \K 4 ) , \K 5 ) , \K e ) , \K 7 ) , \K 8 )} 



(20) 
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which can be explicitly written as 
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-1 
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-1 


-1 


-1 
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-1 


-1 
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-1 
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2 





-1 


-1 
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-1 


-1 


2 
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1 


1 


1 


1 


-1 


2 





-1 


1 


-1 


1 



(21) 



It possesses a number of interesting properties including a computationally useful fact that 
KK) is a diagonal matrix. 



III. CONTINUOUS CENTRAL MOMENTS: EQUILIBRIUM AND SOURCES 

Consider an athermal fluid in motion which is characterized by its local hydrodynamic 
fields at the Cartesian coordinate (x,y), i.e. density p, hydrodynamic velocity it = (u x ,u y ), 
and subjected to a force field F* = (F x , F y ), whose origin could be either internal or external 
to the system. The local Maxwell-Boltzmann distribution, or, simply, the Maxwellian in 
continuous particle velocity space (£ x ,£y) is given by 

2" 



P 



2nd 



exp 



(?-*)■ 



where we choose 



a] = 1/3. 



(22) 



(23) 



Let us now define continuous central moments, i.e. moments displaced by the local 
hydrodynamic velocity, of order (m + n): 



/oo poo 
/ f M (L-u x ) m (t y -u y y 
-oo J —oo 



di x d^ y . 



(24) 



By virtue of the fact that f M being an even function, H^yn ^ when m and n are 
even and H.^ yn = when m or n odd. Here and henceforth, the subscripts x m y n mean 
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xxx ■ ■ ■ m — times and yyy • ■ ■ n — times. Thus, evaluating this quantity in the increasing 
order of moments gives 



r^ 1 = P 
= o 



X 



IF* = 

n y u, 

^xx ~ C sPi 



1L yy 



iF 1 = o 

xy w ' 

fi M = o 

xxy v 

= o 

x yy w ' 

^ = C 4 n 
xxyy ^sr' 

Here, and in the rest of this paper, the use of "hat" over a symbol represents values in the 
space of moments. 

Now, we propose that the continuous distribution function / is modified by the presence 
of a force field as given by the following ansatz: 



Af F = -- K ^ - UJ f M (25) 

P Cs 



It may be noted that He et al. (1998) [38| proposed similar form for the continuous Boltz- 
mann equation to derive source terms for the SRT-LBE. However, it's influence on discrete 
distribution function due to cascaded collision process via the method of central moments 
to establish Galilean invariant solutions is expected to be, in general, be different. Let us 
now define a corresponding continuous central moment of order (m + n) due to change in 
the distribution function as a result of a force field as 



OO /"OO 



r' = / / A/ - U X ) m {iy ~ UyTd^y. (26) 



OO J — oo 



12 



Evaluation of Eq. ( )26|) in the increasing order of moments yields 



1 o — u ? 

f F = F 

f F = F 

f F = 

1 XX u i 

f F = 

1 2/2/ 

f F = 

f F = c 2 F 



xyy 



V = 0. 

xxyy " 



IV. CASCADED LATTICE-BOLTZMANN METHOD WITH FORCING TERMS 



First, let us define a discrete distribution function supported by the discrete particle 

f =|/«> = (/o,/i,/2,...,/8) t . (27) 



velocity set e a : 



Following Geier et al. 32j, we represent collision as a cascaded process in which the effect of 
collision on lower order moments successively influences those of higher order in a cascaded 
manner. In particular, we model the change in discrete distribution due to collision as 

n£ = n£(f,g) = (£•§)„, (28) 

where 

g = \9a) = (do, di,92, • • • , ?s) t (29) 

determines the changes in discrete moment space in a cascaded manner. That is, in general, 

9a = 9a^,9f}), 13 = 0,1,2,..., a- 1. (30) 

The detailed structure of g will be determined later in Sec. EJ 

We define that f a changes due to external force field by the discrete source term S a . 
That is, 

S=\S a ) = (S ,S 1 ,S 2 ,...,S 8 y. (31) 
13 



We suppose that particle populations are continuously affected by this in time as they 
traverse along their characteristics. The precise form of S a is yet unknown and will be 
determined as part of the procedure presented in Sec. [Vj 

With the above definitions, the evolution of f a in the Cascaded-LBM can be written as 

rt+l 

fatf +-?«,t + l) = f a tf, t) + Q c a ^ >t) + J S a & +1 t ag>t+0) de, (32) 

where the fluid dynamical variables are determined by 

8 

P =^fa={fa\p), (33) 
a=0 
8 

pui = ^faCai = (f a \e ai ) ,i <E x,y. (34) 

The last term on the right-hand-side (RHS) of Eq. (132]) represents the cumulative effect of 
forces as particle populations advect along their characteristic directions. Various approaches 
are possible here to numerically represent this integral, with the simplest being an explicit 
rule. However, in general cases where F can have spatial and temporal dependencies, for 
improved accuracy, it becomes imperative to represent it with a higher order scheme. One 
common approach, which is employed here, is to apply a second-order trapezoidal rule, 
which will sample both the temporal end points, (t, t + 1), along the characteristic direction 
a. That is, 

f a {£ + ~i a , t + 1) = / a ( 2^, t) + ft£(-£ >t ) + - [S a ($,t) + S a(-£+l? a ,t+l)] ( 35 ) 

Equation fl35|) is semi- implicit. To remove implicitness along discrete characteristics, we 



apply the following transformation 38|, |42j : 



fa = fa~ \Sa- (36) 

Thus, Eq. (|35|) becomes 

7 a (^ + ~Z m t + 1) = j a tf, t) + + s^. (37) 

Clearly, we need to determine ^ Q S a and J2 a Sa~Za to obtain p and p u , respectively, in 
terms of the transformed variable f a , which will be carried out in the next section. 
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V. CONSTRUCTION OF CASCADED COLLISION OPERATOR AND FORCING 
TERMS 



In order to determine the structure of the cascaded collision operator and the source 
terms in the presence of force fields, we now define the following discrete central moments 
of the distribution functions and source terms, respectively: 



HZ x m y n 



^ ^ faiS'nx ^x) i&ay ^y) ({^qsc ^i) ^ficny ^y) l/a) j ("^) 
a 

0~ x m y n = ^ S a {e ax — U x ) (&ay ~ Uy) = ((&ax ~ U x ) (e ay ~ U y ) \S a ) . (39) 
a 

We also define a discrete central moment in terms of transformed distribution function to 
facilitate subsequent calculations: 



fai^ax ~ Ux) m (e ay - u y ) n = ((e ax - u x ) m (e ay - u y ) n \f a ) . (40) 

a 

Owing to Eq. f !5B|) . it follows that 

Let us also suppose that f a and f a have certain local equilibrium states represented by 
f* q and f Q , respectively, and the corresponding central moments are 

^ ] fa( e ax — u x ) [e ay — Uy) = ((e ax — u x ) (e ay — u y ) \f a q ) , (42) 

a 

J27^(e ax ~ u x ) m (e ay - u y ) n = ((e ax - u x ) m (e ay - u y ) n \J^) . (43) 



Now, we take an important step by equating the discrete central moments for both the 
distribution functions (equilibrium) and source terms, defined above, with the continuous 
central moments derived in Sec. IIHI Thus, we have 

-eg _ frM (AA\ 

l% x m y n — ll x m y n, l^^J 
0~ x m y n = Y x m y n- (^5) 

In other words, the discrete central moments of various orders for both the distribution 
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functions (equilibrium) and source terms, respectively, become 



and 



^eq 

Kq 


= Pi 


K eq 

X 


= o, 


n.. 

y 


= o, 


XX 


= c Ip, 


^yy 


= c 2 s p, 


K'xy 


= o, 


n xxy 


= o, 


K xyy 


= o, 


K 

xxyy 


4 

= c s p, 


?0 


= o, 




— f 


y 


— F 

u 1 


@xx 


= o, 


a yy 


n 

— u, 


@xy 


= o, 


@xxy 


- c 2 F 


®xyy 


= c 2 F 

X) 


@xxyy 


= 0. 



(46) 
(47) 
(48) 
(49) 
(50) 
(51) 
(52) 
(53) 
(54) 



(55) 
(56) 
(57) 
(58) 
(59) 
(60) 
(61) 
(62) 
(63) 

From Eq. (141 p . we get the following transformed central moments, which comprises as one 
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of the main elements for subsequent development and analysis: 



Hn 



^eq 



K 



^eq 

K yy 

^eq 
K xy 



^eq 
^xxy 



^eq 
K xyy 

^eq 
K 

xxyy 



2 

C 2 S P, 



c 2 s P, 



0, 

c 2 
--F 



--F 

_ T. 



C 4 S P- 



(64) 
(65) 

(66) 
(67) 
(68) 
(69) 
(70) 

(71) 
(72) 



To proceed further, we need to obtain the corresponding moments in rest or lattice frame 
of reference, i.e. raw moments. The tool that we employ for this purpose is the binomial 
theorem. The transformation between the central moments and the raw moments for any 
state variable ip supported by discrete particle velocity set can be formally written as 



((e ax -u x r(e ay -u y T\v) = (CC» + (e 

m 
_Z=] 



^C™e™j/(-l) J w^ 
,i=i 



ay 



i c ai I ± ) u x 



m e~ l (-i)X 



i=i 



£W(-i: 

.i=i 



\<P) (73) 



where = p\/{q\{p — <?)!)• In the above, commutation of the inner product of vectors, 
represented using the "bra-ket" operators, with summations and scalar products is assumed. 
Clearly, raw moments of equal or lesser order in combination is equivalent to central moments 
of a given order. 

Application of Eq. (1731) to the forcing terms, i.e., using Eq. (|39|) and Eqs. fl55|) - fl63|) yields 
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analytical expressions in the rest frame of reference: 

(S«\p) = ^S a = 0, (74) 

a: 

{S a \e ax ) = ^ S a c ax = F x , (75) 

a 

(S a |eL> = ^2S a e 2 ax = 2F x u x , (77) 

a 
a 

(<Sq! \&ax&ay) ^ ^ S a C ax &ay F x U y F y U x , (^9) 

a 

(S^e^e^) = Yl S ^x^ y = (l + u l) F v + 2F x u x u y , (80) 

(Sale^e^) = ^ Sae^e^y = ( - + u 2 ] F x + 2F y u y u x , (81) 

(S a \e 2 ax el y ) = S a e 2 ax e 2 ay = Q + 2^J + (j + 2^) F y u y . (82) 

For subsequent procedure, we also need the raw moments of the collision kernel 

£(/C • g) a e™ eS, = ( K ?K*e n ay) 9p- (83) 

Since collisions do not change mass and momenta, which are thus called collisional invariants, 
we can set 

9o = 9i= 92- (84) 

Thus, we effectively need to determine the functional expressions for "gp for — 3,4, ... ,8. 
Owing to the orthogonal property of the eigenvectors of K. by construction, i.e. Eq. f !2Up . we 
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obtain 

X)(^-k)« = Z)(^Ip>^ = °> ( 85 ) 

a p 

22(JC ■ g) a e ax = 22 ( K p\ e ax) dp = 0, (86) 

a /3 

J^(/C ■ g) a e ay = ^ ( K p\ e ay) 9p = 0, (87) 

a p 

]T(/C ■ g) a e 2 ax = Yl ( K P\*L) 9p = 6? 3 + 2? 4 , (88) 
J2(JC ■ g) a e 2 ay = (KpK v ) 9P = 6? 3 - 2? 4 , (89) 
^(AC • g) a e ax e ay = 22 ( K p\ e a X e ay ) gp = 4g 5 , (90) 

a P 

J^(£ ■ S) a e 2 ax e ay = Y ( K p\ e L e ay ) dp = -4? 6 , (91) 

a p 

■ g) a e ax e^ = ^ (i^e^e^) ^ = -Ag 7 , (92) 
^2()C ■ g) a e 2 ax e 2 ay = Y ( K p\ e L e l y ) 9P = 8? 3 + 4# 8 . (93) 

a p 

Now, for computational convenience, the evolution equation, Eq. (|37|) . of the Cascaded- 
LBM with forcing term may be rewritten as 

L(^) = 7atf,t) + n^ >4) + s a{ ^ t) , (94) 
7 Q (^ + ^ a ,t + i) = (95) 

where Eq. AMI) and Eq. (|95|) represent the collision step, augmented by forcing term, and 
streaming step, respectively. Here and henceforth, the symbol "tilde" (~) refers to the 
post-collision state. The hydrodynamic variables can then be obtained as 

8 

p = E7« = <7». ( 96 ) 

8 - 1 - 1 

pui = 22f a e<xi + = (faM + ex ^y ( 97 ) 

a=0 

in view of Eqs. flU}, (JZU), ([75]) and ([76]). 

Now, to obtain the source terms in particle velocity space, we first compute (Kp\S a ), 
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P = 0, 1, 2, . . . , 8. From Eqs. ([H and flUD-flED, we readily get 



m* = 


(K \S a ) 


= o, 


(98) 


ml = 




= F 


(99) 


m s 2 = 


(K 2 \S a ) 


= F 


(100) 


ml = 


(K 3 \S a ) 


= Q(F X U X + FyUy), 


(101) 


fh% = 


(Ki\S a ) 


^{-^X^X -^V^U/l 


(102) 


ml = 


(K 5 \S a ) 




(103) 


rh s 6 = 


(K 6 \S a ) 


= (1 - 3u 2 x )F y - 6F x u x u y , 


(104) 


m s 7 = 


(K 7 \S a ) 


= (1 - 3U 2 y )F x - 6FyUyU X , 


(105) 


m s 8 = 


(K 8 \S a ) 


= 3 [(6uJ - 2)F x . Ml . + (6u* - 2)F y Uy] . 


(106) 



Thus, we can write 

(JC ■ S) a = «#o|Sa) , (Kx\S a ) , (K 2 |S a > , . . . , (K 8 \S a )) 

= (m s , m*, m s 2 ,..., m*) T = |m* ) . (107) 

By virtue of orthogonality of JC, we have KX) = D = 

di ag ((K \K ) , (K X \K X ) , (K 2 \K 2 ) , . . . , (K 8 \K 8 )) = diag(9, 6, 6, 36, 4, 4, 12, 12, 36). In- 
verting Eq. (11071) by making use of the property K~ x = JC* ■ D^ 1 , we get explicit expressions 
for S a in terms of F and it in particle velocity space as 



So = 


- [-ml + m s 8 ] , 


(108) 


Si = 


^ [6m{ — m*. + 9m\ + 6m 7 — 2m*,] , 


(109) 


S 2 = 


[6m* -m s 3 - 9m* + 6m* - 2m*] , 


(110) 


S 3 = 


^ [-6m* - m*. + 9m s 4 - 6m*. - 2m*,] , 


(111) 


S, = 


[-6m* - m* - 9m* - 6m* - 2m*] , 

36 


(112) 


s 5 = 


[6m* + 6m?, + 2m*. + 9m*, - 3ml - 3rh 7 + m s 8 ] , 


(113) 


s 6 = 


- 1 - [-6m* + 6m*, + 2m s 3 - 9rh s 5 - 3rh s 6 + 3m 7 + rh s 8 ] , 
36 


(114) 


s 7 = 


— [-6m* - 6m*, + 2m*, + 9rh s 5 + 3m s 6 + 3rh 7 + rh s 8 ] , 


(115) 


S 8 = 


[6mj - 6rh 2 + 2m*. - 9m*, + 3m s 6 - 3m s 7 + m s 8 ] . 


(116) 
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We now need to find the expressions of {f a \e™ x e^ y ) = Y^ a =o fa e ™x e ay to proceed further. 
In this regard, for convenience, we define the following notation for a compact summation 
operator acting on the transformed distribution function f a : 

a(7 ai +/«,+/«, + • ••)+&(7/j 1 +7 a +7 a + - ••) + •••= + W- ( 117 ) 



Q 



where A = {a±, «2, «3, • • • }, B — {f3\, 02, Ps, • • • },• • • ■ For conserved basis vectors, we have 
them in terms of collisional invariants 

8 



</» = £/« = ft (H8) 

8 - 1 

(f a \e ax ) = E fa e *x = pu x - -F x , (119) 

a=0 

8 - 1 

(faKy) = E f^ay = pUy ~ ~F y , (120) 

a=0 

and, for the non-conserved basis vectors, we have 

(7JO = YJJL = |E ] (121) 



a=0 



A 4 



</ a i4) = E^4 = E ®/«' ( 122 ) 

q=0 \ a ) 

(7je Q;E e a? ,) = E7 Q e CM; e a j / = f E ~ E ) ® 7 a , (I 23 ) 

a=0 \ a a / 

8 / A 6 B 6 \ 

(7 a \e 2 ax e ay ) = E7«^,e ay = I E " E I ® 7«, (124) 

a=0 V a a / 

8 / A 7 B 7 \ 

(7a|e«xO = E7« e «4y = I E " E ) ® 7«, (125) 

8 / A 8 \ 

(7aleL4> = E7a^L4 = E ®7«, (126) 



a=0 
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where 



^3 


= {1,3,5,6,7,. 


3}, 


(127) 


A 4 


= {2,4,5,6,7,. 


s}, 


(128) 


A 5 


= {5,7}, B 5 = 


{6,8}, 


(129) 


A e 


= {5,6},£? 6 = 


{7,8}, 


(130) 


A 7 


= {5,8}, B 7 = 


{6,7}, 


(131) 


A 8 


= {5,6,7,8}. 




(132) 



With the above preliminaries, we are now in a position to determine the structure of the 
cascaded collision operator in the presence of forcing terms. Starting from the lowest order 
non-conservative post-collision central moments, we successively set them equal to their 
corresponding equilibrium states. Once the expressions for "gp is determined, we discard 
this equilibrium assumption and multiply it with a corresponding relaxation parameter to 
allow for a relaxation process during collision [32]. From Eq. (|67|) . which is the lowest non- 
conserved central moment, and applying the binomial theorem (Eq. (1731) ) to transform it to 
the rest frame of reference, we get 

Kl = l/3p = (7jeL> - 2 u x (f a \e ax ) + u\ (f a \p) . (133) 

From Eq. fl94|) and substituting for various expressions involving 
(f a \ e ax)> J2p (Kp\ e ax) 9p and (S a \e™ x ) , where m = 0, 1, 2 from the above, yields 

6^ + 2? 4 =\p- (yyj ® /« + p u * - p * u *- ( i34 ) 

Similarly, from Eq. (I68j) 

Kl = 1 / 3 P = (falely) - 2u y (f a \e ay ) + u\ (f a \p) , (135) 

and using (/Je™ ), Y^p i^/3\ e a y ) 9p an d (^cje™ ), where m = 0, 1, 2 from the above, via the 
binomial theorem gives 

6? 3 - 2g 4 = lp- N~^J ® J a + pu 2 y - F y u y . (136) 
Solving Eq. f |134|) and (I136p for g 3 and g 4 yields 



- 1 2 
93 = i2 3 P 



J2+2j2)®fa + P{< + u 2 y ) - (F x u x + F y u y ) \ , (137) 
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and 

E4, F4, 



where 

C 3 = {1,2,3,4}, (139) 

D 3 = {5,6,7,8}, (140) 

E 3 = {2,4}, (141) 

F 3 = {1,3}. (142) 

Now, we drop the assumption of equilibration considered above applying relaxation param- 
eters, uj 3 and Co> 4 , to Eq. f !137j) and (I138p . respectively, to get 

?3 = ^3^|- fe+ 2 Xj ®7* + lp + p( u l + u2 y)-( F * u z + F y u y)Y ( 143 ) 



and 

94 = u A 



\{(t-t)*f* + p ^ ~ ^ ~ ( F * Ux - F ^ } ' (144) 



Let us now consider the central moment k x in Eq. (169]) . i.e., 



>C = = (f a \{e ax - u x ){e ay -u y )), (145) 



and substituting the expressions for various raw moments, we get 

95 = J |- ~~ Yjj ® 7a + PUxUy - ^(F x u y + F y u x )^j , (146) 

and applying a corresponding relaxation parameter to represent over-relaxation for this 
moment, we obtain, 

95 = ^ j- ~ Y^j ® 7 a + PU x u y - ^(i 5 ^ + i^^)| • (147) 

It is worth noting that due to a slightly different choice of the basis vector K§ for \e ax e ay ) 
from that in |32] , Eq. ( I147p differs from that in |32] by a factor of — 1 apart from the presence 
of forcing terms. 

We now consider the central moment of the next higher order, i.e. 7c^ in Eq. (|70p . 
K xxy = —\F y = (f a \(e ax — u x ) 2 (e ay — u y )) and following the procedure as discussed above, 
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we get 



96 



1 

4 



A & B 6 

E-EHME-E 



it, 



E 



fa 



+2pu x u y + -(1 - u x )F y - F x u x u y > - 2u x g 5 - -u y (3g 3 + g 4 ) 



(148) 



Notice that g$ depends on gp, (3 < 6, which are already post-collision states. So, we relax 
with relaxation parameter uq only those terms that do no contain these terms, leading to 



1 



UJ & - 



A 6 B 6 



A 5 B 5 

E-E)-*. (E-E 



it, 



A 3 ■ 

E 



a a 



fc 



+2pu 2 x u y + ^(1 - ul)F y - F x u x Uy^ - 2u x g 5 - ^- Uy (3g 3 + <7 4 ), 



(149) 



That is, g 6 = g e ({f a } , p, it, ~P , g 3 , 94, 9b, We)- 

Considering next, = —\F X = (f a \(e ax - M :c )(e Qy - u y ) 2 ) from Eq. (|7TJ) and following 
calculations to transform all the quantities to raw moments, we get 

A 5 B 5 
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1 

4 



A 7 B 7 

E-E 



-2u, 



E-E 



A 4 ■ 



fc 



-2pu x u 2 y + ^(1 - u y )F x - FyU y u x j - 2w y <?5 - ^m x (3^ 3 - g A ), 



(150) 



Again, notice that g-j depends on gp, (3 < 6, which are already post-collision states. So, 
applying the respective relaxation parameter u-j to terms that do no contain them, yields 



1 

97 = u 7 - 



A 7 B 7 



E-E E-E -«*E 



/a 



+2pu x u 2 y + ^(1 - mJ)F z - FyU y u x X - 2u y g 5 - ^u x (3g 3 - </ 4 ), (151) 

Thus, #7 = <77({/ a } , p, ~it, ^,?3,?4,^5, w 7 ). In other words, depends on only the lower 
order moments and not on other components of the same order. 

Finally, we consider the central moment of the highest order defined by the discrete 

P = (f a \( e ax ~ u x ) 2 (e ay - u y ) 2 ), and apply the 



particle velocity set (Eq. (172]) ). k xx 



yy 9 
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procedure as discussed above to transform everything in terms of raw moments to obtain 

A A A 3 



1 

4 

1 r 



A s / A 7 B 7 

E- 2 «dE-E 



2u, 



A 6 B 6 

E-E 



u: 



As 



-Br, 



E-E 



/ Q + -p + 3pu*uJ - (F x u x Uy + F y u y u 2 x ) 



2?3 



?<0 - ^m^(3?3 + 9i) ~ 4u x ii y g 5 - 2u y % - 2u x g 7 , 



(152) 



Clearly, g% depends on "gp, j3 < 7, which are already post-collision states and thus, we relax 
with the parameter u$ those terms that do not contain them to finally yield 



W8 I 



4:U x 1ly 



'As / A 7 B 7 

E- 2 ^ E-E 



A 6 



Be 



A 4 



2m,. 



Br, 



E-E 



E-EME+sE+ 

-2? 3 



1 fa + \p + 3p?4?4 - + 



94) 



1 2 
2^ 



+ 9a) ~ ^u x u y g b - 2u y g 6 - 2u x g 7 , 



(153) 



In order words, g 8 = g 8 ({f a } , p, ~it, ?3, 94, ?5, 9e, 97, w 8 ). It may be noted that because of 
a slightly different choice of the basis vector K 5 , the prefactors for g 5 in Eqs. fll49p - fll53p 
differ from that in [32J by —1. Unfortunately, in the seminal work 32], there are some 
typographical errors in Eqs. (20)- (24) of that paper jsjj - in particular, some of the signs 
in the last lines of its Eq. (20)- (23), and the expression in the last line of its Eq. (24) are 
incorrect. 

Thus, the general structure of cascaded collision operator for non-conserved moments 
may be written as 

g a = u a it) *M({fp}) + H 2 (j>, it) o N(f)} + C(g,), 



(154) 



where a = 3, . . . , 8, (3 = 0, 1, 2, . . . , 8 and 7 = 0, 1, 2, . . . , a — 1, and M, N, Hi, and H 2 
represent certain functions, and * and o represent certain operators. On the other hand, 
in particular, the term C(g 7 ) contains the dependence of g a on its corresponding lower 
order moments leading to a cascaded structure. In other words, cascaded collision operator 
markedly distinguishes from the SRT and MRT collision operators in that the former is non- 
commutative. The above derivation involved the choice of a particular form of the central 
moments of the sources. In the next section (Sec. IVII) . it will be shown how a different choice 
could provide a better representation of its effect on higher order moments. 
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VI. DE-ALIASING HIGHER ORDER CENTRAL SOURCE MOMENTS 



Due to the specific formulation of the forcing term employed in Eq. (|25|) . its corresponding 
higher order central moments also have non-zero contributions, even when the fluid is at 
rest and a homogeneous force is considered. Since they only occur at third and higher 
order moments, they do not affect consistency to the Navier-Stokes equations, which emerge 
at the second-order level (see Appendix [A]) . However, to be conceptually consistent, it is 
desirable to avoid this effect. Thus, as a limiting case, we now maintain the effect of the 
force field only on the components of the first-order central source moments, and de-alias 
all the corresponding higher (odd) order central moments, by setting them to zero. That is, 



1 x m y n 



F x , m — 1 , n — 

F y , m = 0,n = 1 (155) 
0, m + n > 1. 



In effect, the transformed equilibrium central moments K xmy n used in the construction of the 
collision operator are modified. Specifically, the third-order transformed equilibrium central 
moments, Eqs. (170]) and ( I7TT) now reduce to 

K xxy = K xyy = 0' (156) 

while all the other components are the same as before. Moreover, such de-aliasing also 
modifies the raw moments of the forcing terms at higher orders. In particular, Eqs. ( 1801) - 
( 1521) now become 

(Sa\& ax e ay) = ^ ] S a e ax e ay = F y U x + 2F x U x U y , (157) 

a 

(S a \e ax e 2 ay ) = ^ S a e ax e 2 ay = F x u 2 y + 2F y u y u x , (158) 

a 

(S a \e 2 ax el y ) = s a e 2 QX e 2 ay = 2F x u x u 2 y + 2F y u y u 2 x . (159) 

a 

while the lower order moments remain unaltered. Notice that terms such as 1/3F X and l/3F y 
do not anymore appear in the third-order source moments, while 2/3F x u x and 2/3F y u y are 
eliminated from the fourth-order source moments as a result of the use of de-aliased central 
source moments (Eq. (j!55p ). Hence, when the fluid is rest, the force fields do not influence 
the third and higher order raw source moments, which is physically consistent. 
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The computation of the source terms in velocity space S a using Eqs. fll08l) - (lll6l) . which 
involve rhi, are also naturally influenced by the above changes. In this regard, while m^, for 
(3 = 0, 1, 2, . . . , 5 remain unmodified, the higher order moments for /3 = 6, 7, 8 are altered. 
The expressions for these latter quantities now become 



m. 



(K 6 \S a ) = (2 



3u 2 x )F y 



m s 7 = (K 7 \S a ) 
ml = (K 8 \S a ) 



(2 - 3u 2 y )F x - 6F y u y u x , 

6 [(3uJ - 2)F x u x + (3u 2 x - 2)F y u y ] 



(160) 
(161) 
(162) 



The cascaded collision operator can now be constructed using the procedure presented in 
Sec. |Vj The use of modified source moments do not alter the collision kernel corresponding 
to gp, where (3 = 0,1,2,. ..,5 and (3 = 8. They are the same as those presented in Sec. [V] 
On the other hand, the third-order collision kernel contributions are modified, which are 
now summarized as follows: 



96 



W6 4 




Br, 



- 2w, 



E-E -«.£ 



fa 



+2pu 2 x u y - -u\F y - F x u x u y \ - lurf}*, - ^(3^ + g A ), 



(163) 



and 



97 



uj 7 - 




Br, 



2u,. 



E-E 



A 4 

u x J2 



+2pu x u 2 y - -u 2 y F x - F y u y u x \ - 2u y g 5 - -u x 



3-94)- 



(164) 



Again, evidently, when the fluid is at rest, the force fields do not have direct influence on 
^6 and 'g-j. Thus, the above formulation eliminates spurious effects resulting from forcing 
due to the finiteness of the lattice set for higher order moments, similar to that by Guo et 
al. 4l| for other LBM approaches. Indeed, a Chapman- Enskog multiscale moment expansion 



analysis carried out in Appendix |A] will establish the consistency of this special formulation 
of the central moments based LBM to the desired macroscopic fluid flow equations. The 
shear and bulk kinematic viscosities is found to be dependent on the relaxation parameters 
w 3 = u x and u>4 = 00$ = u u , respectively. In particular, the shear viscosity satisfies v = 
c s — l) ■ The rest of the relaxation parameters in this MRT cascaded formulation can 
be tuned to maintain numerical stability. One particular choice suggested by Geier is to 
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equilibrate higher order, in particular, the third-order moments, 1 |35|. 

Other possible choices could be also considered that involve over-relaxation of these moments 
at certain carefully selected relaxation rates so as to control numerical dissipation while 
maintaining computational stability. On the other hand, as shown in Appendix [Bj when the 
central moments based LBM as derived in this work is executed as a MRT cascaded process 
it implies generalization of both equilibrium and sources in the lattice frame reference which 
also depend on the ratio of various relaxation times. However, it does not affect the overall 
consistency of the approach to the macroscopic equations as it influences only higher order 
contributions. The discussions so far considered the cascaded collision operator to be explicit 
in time. Appendix O presents with the possibility of introducing time-implicitness in the 
cascaded collision operator. 



VII. COMPUTATIONAL PROCEDURE 

The main element of the computational procedure consists of performing the cascaded 
collision, including the forcing terms, i.e. Eq. ( 193)) along with Eq. ( |28|) . which can be 



expanded as follows: 

70 = 7o + [?o - 4(? 3 - 9s)] + S , (165) 

71 = 7 1 + [?0+?l-#3+?4 + 2(?7-?8)] + S 1 , (166) 

72 = 72 + [do + 92 -93 -94 + 2(? 6 - ds)} + S 2 , (167) 

73 = 73 + [do - 91 - 93 + 94 - 2(<? 7 + 0b)] + S 3 , (168) 

7 4 = 7 4 + [do - 92 - 93 - 94- 2(# 6 + g&)] + 5*4, (169) 

7 5 = 7 5 + [do + 91+92 + 2<? 3 + 95-96-97 + 9s] + S B , (170) 

7 6 = 7 6 + [do -91 + 92 + 2? 3 -95-96 + 97 + 9s] + Sq, (171) 

77 = 77 + [do -91-92 + 2<?3 + 95 + 96 + 97 + 9s] + S r , (172) 
7s = 7 8 + [do + 9i-92 + 2? 3 - 95 + 96 - 97 + 9s] + S 8 . (173) 



Here, the terms "gp can be obtained in a sequential manner, i.e. evolving towards higher 
moment orders from Eqs. ( 1135]) . f lT44j) . ffMTj) . ffT49|) . ffToTj) . and ffT53l) . It consists of terms 
that involve summation of f a over various subsets of the particle velocity set. The source 
terms Sp are computed from Eqs. (1108)) - fjl 161) . Once the post-collision values, i.e. f a are 
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known, the streaming step can be performed in the usual manner to obtain the updated 
value of f a (Eq. fl95l) ). Subsequently, the hydrodynamic fields, viz., the local fluid density 
and velocity can be computed from Eqs. ( 1961) and ( 1971) . respectively. Depending on the 
specific choice of the ansatz for the central source moments, appropriate expressions for gp 
and rhp need to be used (see Sees. IVl and IVTl) . In the above procedure, careful optimization 
needs to be carried out to reduce the number of floating-point operations. 

VIII. COMPUTATIONAL EXPERIMENTS 

In order to validate the numerical accuracy of the new computational approach presented 
in this work, we performed simulations for canonical fluid flow problems subjected to different 
types of forces, where analytical solutions are available. We will now present results obtained 
by employing the Cascaded-LBM with de-aliased higher order source central moments (as 
discussed in Sec. I VI I) , which will be compared with corresponding analytical solutions. The 
first problem considered is the flow between parallel plates subjected to a constant body 
force. We considered 3 x 51 lattice nodes to resolve the computational domain, where 
periodic boundary conditions are imposed in the flow direction and the no slip boundary 
condition at the walls is represented by means of the standard link bounce back technique. 
The relaxation parameters are given such that = u§ = 1.754, while the remaining ones 
are set to unity and the computations are performed for different values of the component 
of the body force in the flow direction, i.e. F x with F y = 0. Figure [2] shows a comparison of 
the computed velocity profiles with the standard analytical solution (Poiseuille's parabolic 
profile, with the maximum velocity uq = F X L 2 /(2z/), where L is the half-width between the 
plates and v is the fluid's kinematic viscosity) for different values of F x . Excellent agreement 
is seen. In order to quantify the difference between the computed and analytical solution, 
the relative global error given in terms of the Euclidean (second) norm is presented in Table 
I. Thus, for the above given set of parameters and resolution, it is O(10 -4 ). 

The second problem considered involves a spatially varying body force. One classical 
problem in this regard is the Hartmann flow, i.e. flow between parallel plates subjected to 
a magnetic field B y = Bq imposed in the perpendicular direction to the fluid motion. If Ff, 
is the driving force of the fluid due to imposed pressure gradient and Ha is the Hartmann 
number that characterizes the ratio of force due to magnetic field and the viscous force, then 
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15 20 25 



FIG. 2: Flow between parallel plates with constant body force: Comparison of velocity profiles 
computed by Cascaded-LBM with forcing term (symbols) with analytical solution (lines) for dif- 
ferent values of the body force F x . 



Magnitude of body force (F x ) 


Relative global error ( \\5u 2) 


1 x 10~ 6 


3.999 x 10" 4 


3 x 1(T 6 


3.895 x 10 -4 


5 x 1(T 6 


3.837 x 1(T 4 


7 x 1(T 6 


3.839 x 1(T 4 



TABLE I: Relative global error for the Poiseuille flow problem. 1 1 <5xi 1 1 2 = SJK^c.i — 
U a,i)\\2/Yli H u o,i||2) where u Cj i and u a ^ are computed and analytical solutions, respectively, and 
the summation is over the entire domain. 



the induced magnetic field in the flow direction B x is given by B x 



F b L 
B 



sink (Ha-^) 



stn/i(Ha) L 

where the coordinate distance y is measured from a position equidistant between the plates. 
The interaction of the flow field with the magnetic field results in a variable retarding 
force F mx = B y <^ and F - " dB 
is F x = F b + F mx and F y 

the same values of the relaxation parameters as before, with F^ = 5 x 1CT 6 and B 



= —B x ^-, and, in turn, the net force acting on the fluid 
F my . We considered the same number of lattice nodes and 
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FIG. 3: Flow between parallel plates with a spatially varying body force: Comparison of velocity 
profiles computed by Cascaded-LBM with forcing term (symbols) with analytical solution (lines) 
for prescribed Lorentz force at different Hartmann numbers. 



8 x 10 3 and varied the values of Ha. The analytical solution for this problem is it a 



Zgy/IcOthQto) 



cos/i(Ha^- J 
cosh(Ra) 



where the magnetic resistivity rj is related to Ha through 



T) = H ^ ■ The computed velocity profiles are compared with the analytical solution for 
different values of Ha in Fig. |3] As expected, the velocity profiles become more flattened 
with increasing values of Ha, while the case with Ha = reduces to the earlier problem. 
The computed velocity profiles are found to agree very well with the analytical results. The 
relative global errors for this problem are presented in Table II. It can be seen that they are 
dependent on the value of Ha when the same grid resolution is used for different cases. In 
particular, the relative error increases as the value of Ha is increased for the same resolution. 
This can be explained as follows. This flow problem is characterized by the presence of 
boundary layers - the Hartmann layers - whose thickness is inversely proportional to -y/Ha. 
That is, the Hartmann layer becomes thinner as the value of Ha is increased. Thus, resolution 
of this boundary layer would require increasingly more number nodes that are clustered near 
walls as Ha is increased to maintain the same accuracy. Otherwise, when the same number of 
grid nodes that are uniformly distributed is employed, the relatively error norm is expected 
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to increase with Ha. Indeed, local grid refinement employing a suitable boundary layer 
transformation can maintain similar accuracy for different Ha as was done with other LBM 
formulations recently Q|. Extension of the local grid refinement approaches for the central 
moment based LBM to resolve boundary layers and sharp gradients in solutions are subjects 
of future studies. 



Hartmann number (Ha) 


Relative global error ( \\Su 2) 


0.0 


3.837 x 10" 4 


3.0 


2.140 x 10~ 3 


5.0 


5.967 x 10~ 3 


7.0 


1.091 x 10~ 2 



TABLE II: Relative global error for the Hartmann flow problem. ||<5u||2 = SJK^c.i — 
u a,i)\\2/ J2i \\ u a,i\\2, where ti C) j and u a ^ are computed and analytical solutions, respectively, and 
the summation is over the entire domain. 



The last problem that we considered involves a temporally varying body force. An im- 
portant canonical problem in this regard is the flow between two parallel plates driven by 
a force sinusoidally varying in time. That is, we considered F x = Fi>cos(ut), where is 
the peak value of the applied force, while u p = 2n /T is the angular frequency where T is 
the time period. This problem is characterized by Wo = J^-L, a dimensionless number 
arising from its original analysis by Womersley. The analytical velocity profile for this flow 
is u x = 7Z. jl - ^ffi I e** 1 , where 7 = ^/-iWo 2 . We considered F b = 1 x 1(T 5 
and Wo = 12.71, while maintaining the number of lattice nodes and the values of the re- 
laxation parameters to be same as in the first problem. Figure H] shows a comparison of the 
computed velocity profiles with analytical solution for different instants within the duration 
of the time period T of the cycle. Evidently, the new computational approach is able to 
reproduce the complex flow features for this problem involving the presence of Stokes layer 
very well. Table III presents the relative global errors at different instants within the time 
period T, corresponding to those in Fig. HI The relatively differences between computed 
and analytical solutions vary between different time instants. On the other hand, they are 
identical for instants shifted by the half time period implying that the computations are 
able to reproduce temporal variations without any time lag as compared with analytical 
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0.02 




FIG. 4: Flow between parallel plates with a temporally varying body force: Comparison of velocity 
profiles computed by Cascaded-LBM with forcing term (symbols) with analytical solution (lines) 
at different instants within a time period T. 

solutions. 

It may be noted that for all the three benchmark problems presented above, essentially 
same numerical results are obtained when the de-aliasing in the forcing is turned off, i.e. 
expressions presented in Sec. |V] is used. This is because both forms differ only in third and 
higher orders, while they are both consistent at the second order level with the Navier-Stokes 
equations, from which the analytical solutions are derived. It would be interesting to carry 
out detailed numerical error analysis as well as stability analysis of the central moment based 
LBM for different grid resolutions and characteristic parameters, and for various canonical 
flow problems in future investigations. 



IX. SUMMARY AND CONCLUSIONS 

In this paper, we discussed a systematic procedure for the derivation of forcing terms 
based on the central moments in the Cascaded-LBM. The main elements involved in this 
regard are the binomial theorem that relates the central moments and raw moments of vari- 
ous orders and the associated orthogonal properties. The discrete source terms are obtained 
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Time instant (t) 


Relative global 


error ( \\Su 2) 





4.195 


x 


10~ 3 


0.05T 


1.701 


x 


10~ 3 


0.10T 


1.060 


x 


10~ 3 


0.15T 


7.548 


x 


10- 4 


0.20T 


5.906 


x 


10- 4 


0.40T 


1.842 


x 


10~ 3 


0.45T 


4.611 


x 


10- 4 


0.50T 


4.195 


x 


10~ 3 


0.55T 


1.701 


x 


10- 3 


0.60T 


1.060 


x 


10~ 3 


0.65T 


7.548 


X 


10- 4 


0.70T 


5.906 


X 


10- 3 


0.90T 


1.842 


X 


10~ 3 


0.95T 


4.611 


X 


10- 3 



TABLE III: Relative global error for the Womersley flow problem. ||<5u||2(t) = Yli\\{ u c,i(t) ~ 
u a,i(t))\\2/ ^2i\\ua,i(t)\\2, where u c> i(t) and u a ^(t) are computed and analytical solutions, respec- 
tively, at instant t within a time period T and the summation is over the entire domain. 

by matching with the corresponding continuous central moment of a given order. For the 
latter, we consider an ansatz based on the local Maxwell distribution. Its variant involv- 
ing a de-aliased higher order central source moments, which recovers physically consistent 
higher order effects when the fluid is at rest, is also derived. Effectively explicit and tem- 
porally second-order forms of forcing terms are obtained through a transformation of the 
distribution function, which contributes to the cascaded collision. When the values of the 
free parameters in the continuous equilibrium (Maxwell) distribution, i.e. speed of sound 
and those in the orthogonalization process of the moment basis from the discrete velocity 
set are chosen, they completely determine the various coefficients of both the cascaded col- 
lision operator and the source terms. The equilibrium distribution and the source terms 
in velocity space are proper polynomials and contain higher order terms. By construction, 
the source terms are Galilean invariant. It is found that both the equilibrium and source 
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terms generalize when the cascaded formulation is represented as a relaxation process in the 
lattice frame of reference. While the Cascaded-LBM with forcing terms is based on a frame 
invariant kinetic theory, its consistency to the Navier-Stokes equations is shown by means 
of a Chapman- Enskog moment expansion analysis. It is found that the new approach repro- 
duces analytical solutions for canonical problems that involve either constant or spatially 
or temporally varying forces with excellent quantitative accuracy. The approach presented 
in this paper can be extended to other types of lattices such as the D3Q27 model in three 
dimensions 44 1. 



Appendix A: Chapman-Enskog Multiscale Analysis 

In this section, let us perform a Chapman-Enskog analysis of the central moment for- 
mulation of the LBM using the consistent forcing terms derived in Sec. I VI I For ease of 
presentation and analysis, we will make a particular assumption regarding the collision op- 
erator in this section. It will then be pointed out in the next section that relaxing such 
assumption amounting to the use of fully coherent cascaded collision kernel does not affect 
the consistency analysis presented here. First, some preliminaries are provided. In partic- 
ular, we define a transformation matrix corresponding to the following "nominal" moment 
basis on which the analysis is performed: 

T ~ \\P) ' \ e ax) ) \ e ay) j \ e ax e ay) > \ e ax ~ e ay) ) \ e ax e ay) j \ e ax e ay) > \ e ax e a y) 5 | e ax e aj/)] ) (Al) 

It is convenient to carry out the multiscale expansion in terms of various raw moments. 
Thus, we also define the following raw moments, where the superscript "prime" symbol is 
used here and henceforth to designate that the moment is of raw type: 

(A2) 
(A3) 
(A4) 
(A5) 
(A6) 



, / 

K x m y n 


= V f e m e"' = 

a 


\ e ax e ay\fa) 5 


/ 

<J x m y n 


- 9 P m P n 
/ j u a^ax^ay 

a 


( p m n 1 c \ 
Y^ax^ay \ LJ ctl > 




— feq m n _ 
/ j J a °ax ay 

a 


/ m n I f eq\ 
\ ax ay 1 J a 1 


r / 

t\i,j.m yTl 


= V7 e m e" = 

/ j J a^ax^ay 
a 


(e m e n if ) 

\^ax^ay\J al 1 




— \^ ~f eq p m n _ 
/ j J a ^ax ay 


(e m e n \T q ) 

Votx ay \J a 1 
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It follows that K 

We now re-write various different central moments in terms of their corresponding raw 
moments by applying the binomial theorem. First, the non-conserved part of the central 
moments can be written as functions of various raw moments as follows: 

K X x = k xx - pu 2 x + F x u x , (A7) 



Kyy 



K yy ~ P u v + F v u yi ( A8 ) 



1 



K xy = K xy - pu x Uy + -(F x u y + F y u x ), (A9) 

K xxy = K xxy — 2u x K xy — u y K xx + 2pu x u y — -F y u x — F x u x u y , (A10) 

^xyy = ^xyy ^ U y K xy Ux ^yy ^P U x U y "n^x^y Fy u y u xi (AH) 

K X xyy = K X xyy ~~ ^ U x K X yy ~ ^^y^xxy U x K yy U y K xx ^ u x u y^ X y 

-3pu 2 x u 2 + F x u x Uy + FyU y ul. (A12) 
The raw moments of the equilibrium distribution and source terms of various order are: 

K q ' = P, (A13) 

K q ' = Pu x , (A14) 

K y q ' = pu y , (A15) 

Kt = \p + Pul (A16) 

K e y q y = ^P + PU 2 y , (Al7) 

K q y = PUxUy, (A18) 

Kiy = ^PUy + pU 2 x U y , (A19) 

K q w = ^pu x + pu x u 2 y , (A20) 

Kit,,. = \p+ \p(< + U 2 y) + Pulu'y, (A21) 



"xxyy 
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and 





= o, 


(A22) 


r / 

a x 


— f 


(A23) 


V 


= F 

y 1 


(A24) 


I 

XX 




(A25) 


J 

(J 

yy 


— ?/ 


(A 26) 


/ 

a xy 


^x^y ~t~ ^y^x-i 


(A27) 


/ 

® xxy 


^y^x ^^x^x^yi 


(A28) 


® xyy 




(A29) 


® ' xxyy 


= 2F x u x u 2 y + 2F y u y u 2 x , 


(A30) 



respectively. 

In the above notation, the cascaded collision kernel may be more compactly written as 

93 = Y2 1 3 P + p ^ Ux + u y> ~ ( Kxx + K w> ~ t° xx + ° vy) J ' ^ ' 

94 = y |p(«a - U l) ~ ftxx - %y) ~ ~ ' ( A32 ) 

95 = y \^>UxUy ~ \ y ~ Tf'xy} > ( A33 ) 

?6 = y 2 P^X + « xxy - 2^7^ - u y K xx - -a xxy j- - -w y (3^3 + <? 4 ) - 2m^ 5 , (A34) 

97 = ^ \^P U * U l + "*W ~ 2 V^/ - - ^wj ~ ^( 3 #3 - &) - 1Uy%, (A35) 

- w 8 J 1 2 2 ["-' ~' ~' 2- 2- 

= T I gP + 3pu^ - ^ - 2^ - 2u y n xxy + ^ + u y n xx 

+Au x u y K xy - \f' xxyy j - 2g 3 - ^mJ(3^ 3 + - ^u 2 x {3g 3 - g 4 ) 

-iu x u y g 5 - 2u y g 6 - 2u x g 7 . (A36) 

Instead of considering the above collision operator, for now, in what follows, let us special- 
ize the collision term. In this regard, we first re- write the cascaded collision step, Eq. fl94|) . 
using Eq. (I28I) as 

(JC-g) a = (J Q -J a ) + S aj (A37) 
and reduce it by applying the central moment operator ((e ax — u x ) m (e ay — u y ) n \-) on both 
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of its sides. Thus, we get 



^ (( e ax — U x ) m (e Q y — Uy) n \Kp) gp — (fi^m^n — K x m y n) + (J 



x" l y" 



(A38) 



Let us now consider a specific case when the post-collision state is in "equilibrium state". 
In this case, we set 

K x myn = H x myn, <T x m y n = =^ Q ft = CJ Q (A39) 

so that gp takes certain specific values, gt. 

Thus the specialized non-conserved collision kernel can be obtained by expanding the 
LHS of Eq. (1A38P and using Eq. f]A39|) for m + n > 2, which can be written in matrix form 

as 



7 







- ^eq 

K xx ~ 


^xx 


9l 




^eq 

K yy ~ 


Kyy 


9t 




am 

K xy - 


K'xy 


SS 




aeq 
K xxy 


^xxy 






^eq 
xyy 


- Kxyy 


9*8. 




^eq 
n xxyy 


K'xxyy 



(A40) 



where T = J-(lt,t) is a 6 x 6 local frame transformation matrix that depends on the local 
fluid velocity and is given by 



6 
6 


—6u y 
— 6m t 



2 

-2 


-2Uy 

2u v 






4 

—8u x 















-4 
0-4 



(A41) 



[8 + 6(u 2 x + u 2 y )) -2(ul - u 2 y ) 16u x u y 8u y 8u x 4 

__ n 

It may be noted that Eq. ( 1A41f) has entries similar to that given in Ref. [36], except for the 
change in signs in the third column resulting from the specific choice made for constructing 
\K§) in the orthogonalization (Gram-Schmidt) procedure. Now substituting for the expres- 
sions in the RHS of Eq. (IA40I) and inverting it, we get gt in terms of the raw moments, 
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hydro dynamic fields and force fields. It may be written as 



n* 




1 n 4- 1 
18" ' 12 


1 

P\ U x + U y) ~ u\ K xx + K yy) ~ 


^ ( TP n j. Pii 1 
12 \ i u i ~r~ - 1 y^y) 






1 ^„,2 


/ / 


TP a \ 

JPyUy) 


9t 






1 ' i 

— J^xy ~ gi-FxUy + F y U x ) 




Fe 




— J^pUy 


1 1 1 

— JpU x Uy + -l& XX y + jF x U X Uy 


+ gFyU x 


9*7 




~J2 PU X 


- Jpw*w5 + + J F yUyUx 


+ ±F u 2 








hp{ul + u 2 y ) + \pulu 2 y + \{-K, 


+ ^yy) — l^xxyy + <?£2OT _ 



(A42) 



where g xa:ra = + i^Wj/) — \ (F x u x u 2 + F y u y u 2 ). An alternative and a somewhat direct 

procedure to obtain Cg is to invoke the orthogonal properties of the basis vectors \Kp). 
Accordingly, we can write 



9a 



<JZ-fa\Kp) _ (f e a q -fa- l 2 S a \K p ) 



3, 4,5,. ..,8, 



(A43) 



(Kp\K p ) {K P \K P ) 
which gives expressions identical to that given in Eq. ()A42j) . 

Equivalently, for the special case noted above (Eq. (1A39|) ). the collision operator, 
Eq. ( 1 A3 71) . can also be written as K ■ g* = € q — f = f eg — f — |S, which can be inverted to 
yield 

g* = /C- 1 ^-f-is^), (A44) 

where as before the boldface symbols represent the column vectors. Now, we propose to 
"over-relax" the above special system by means of multiple relaxation times (MRT) as a 
representation of collision process. That is, we set 



g = Aj 



(A45) 



where A is a relaxation time matrix. Hence, combining Eqs. (1A44I) and (1A45I) . we can write 
the post-collision state in this MRT formulation as 



f = f + K • e + S = f + /CAg* + S 



1 



Let, 



Hence, 



f + /CA/CT 1 ( f e9 -f - -S ) + S 



A* = /CA/C" 1 . 



f = f + A* (f eq - f ) + ( X - -A* ) S 



(A46) 

(A47) 
(A48) 
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where X is the identity matrix. 

We now define raw moments of distribution functions (including the transformed one), 
equilibrium and sources for convenience as 

?=Tf, f = Tf, f e <? = rf e9 , S = TS, (A49) 

where (■) represents column vectors in (raw) moment space and the transformation matrix 
T is given in Eq. (lAli . That is, 



— /— — — — \ f / *<j 

f ~~ ( f Oi f li f 2i ■■■if?,) ~ ( K 0i ^xi ^yi K xx ^yyi K xx ~ K yyi ^xyi ^xxyi ^xyyi K i 
f ( foi fit f2j ■ ■ ■ i fs ) ( ^0' K xi Kyi K xx Kyyi ^xx K yy ^xy ^xxyi ^xyyi K xxyy ) i 



f eg _ ( Teq 7eq Teq 7bq\ / -eg' - eg ' -eg' -eg' , -eg' £eg' _ -eg' -eg' -eg' -eg' -eg' \ t 

1 " U0 iJl wJ r--)J8 I I "'O J "'x ' "fy ' ™xx * ^yy J ^xx ^yy i ^xy i ^xxyi ™xyyi ^xxyy I i 

S = ( 5*0, S*i, 5*2, . . . , S*8 J = ( (T , cr x , <r y , a xx + a yy , a xx — <r yy , a xy , cr xxy , cr xyy , <r xxyy 



Finally, using Eq. ( 1A49|) , we can rewrite the expressions for the collision and source terms 
in Eq. (1A48I) in terms of (raw) moment space. That is, 



f = f + r 



-i 



-A (f-f'M • • ( J-^A IS 



(A50) 



where A is a diagonal collision matrix given by 

A = TA*T _1 = diag(0, 0, 0, u 3 , u}±, u 5 , u e , u 7 , (A-51) 

It may be noted that from Eq. (1A49I) . we can obtain the discrete equilibrium distribution 
functions and source terms in velocity space by means of the inverse transformation. That 
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is, f eq = 


r -i f e ?) s = r 


Jo — 


4 

9^ 


2 r 2 

3 v x 


f? = 


1 

9^ 


1 

-pu x + - 

3 : 


req _ 
J 2 — 


1 
9 


l : 

nP U V + - 
3 


pea 


1 

9 


1 

oPUx + - 

3 


req 

J A = 


1 

9 


l : 

nP U V + I 

3 


r('(l 

J5 — 


1 

— P - 

36^ 


H y 2 Pu x - 


J 6 


1 


1 

- T2 PU X 


req _ 

h — 




- T2 PU X - 


req _ 

Js — 


i i 

— P H pu x - 


and 







+ pU 2 x U 2 y , 
ll ll 

2 P U l ~ qP( U 1 + U 2 y) - ~PU X U 2 y ~ -pU 2 x U 2 y , 

2 rJj y~ \ p< ^ + M S) _ - \pu 2 x u 2 y , 

\pu 2 x - ^p(w* + u 2 y ) + ~pu x uj - "P^, 

\p u I - \p( u l + M D + - 

l2 pM?/ + l2 p ^ + u2y) + \ pUxUy + \ pu ^ u y + 4^^ + \p u I u % 
j^puy + + u2 y> ~ \ pUxU y + \f m l u v - \p u * u I + 

1 1 , o . 1 1 o 1 2 1 

pU X U y + ~pU x U y , 



12 y 12 y 4'----y 4' ±r-~*- y ^ 

1 1 , , 1 1 o 1 , 1 

12^+12 



P u y + J^P( U 1 + u l) ~ -{pu x u y - -pu\u y + -pu x u 2 y + -pu 2 x u\, 



So 


^F x 1l x ^^ytky H~ ^FxUxliy -\- ^FylXyU^ 












Si 


-\-—F x -\- F x u x ~F x Uy FyiiyUx F x u x Uy ~ 


- FyUyU x , 










s 2 


1 1 2 2 

~^~~^2~^y FyUy ^FyUj. F X U X Uy F X U X Uy ~ 


- F y UyU x , 










s 3 


1 1 2 2 

= --F x + + -F x u y + - F x u x u y - 


' FyllyU x , 










s, 


1 1 2 2 

2 «' v^y 2 ^^^^y F x u x iiy — 


Fy^ j yU j X '> 




2 F X U X Uy 




-FyUyU x , 


s 5 


1 1 1 2 1^2 l n 

= +^-^% + ^j/Mx + ^F x u w + —F y u x + -F x u 


X Uy ~\~ — FyUyUx 


+ 




s. 


1 1 1^ 2 1^2 l n 

= -^ F xUy - -F y u x - -F x u y + -F y u x + 


X Uy —FyUyU X 


+ 


^F X U X U 2 y 




-FyUyU x , 


s 7 


1 1 1 2 1 2 1 

4 x y ^ y x ^ X y 4 2 


1 

X Uy —J7yUyU X 


+ 


— F X U X Uy 




-FyUyU x , 


s 8 


1 1 1 2 l p 2 l p 

— —-^ x Uy — —r y U x + —r x U y — —r y U x — — r x U 


X Tly H - —FyUyU X 


+ 


^F X U X Ul 




~FyUyU x . 



Thus, the discrete equilibrium distribution and forcing terms in velocity space resulting from 
corresponding imposed central moments are proper polynomials containing higher order 
terms as compared to the standard LBM. The specific functional expressions for f^ g and S a 
depend on the choice made for the "nominal moment basis" (Eq. (1Alj) ) from which they are 
derived. 
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We are now in a position to perform a Chap man- Enskog multiscale expansion. First, 
expand the raw moments f (untransformed ones, i.e. without "overbar", for simplicity) and 
the time derivative in terms of a small bookkeeping perturbation parameter e (which will be 
set to 1 at the end of the analysis) [421 ]: 



f = ^ e n f (n) , 

n=0 

oo 

dt = 



(A52) 
(A53) 



n=0 



We use a Taylor expansion for the representation of the streaming operator, which is carried 
out in its natural velocity space: 



n n 

f + ^ a e, t + e) = V -Ad t + -f a ■ ^)f (^, t). 

z — ' n\ 



(A54) 



n=0 



Substituting all the above three expansions in the LBE, with Eq. ( 1A50|) representing the 
post-collision, and equating terms of the same order of successive powers of e after making 
use of Eq. (1A49j) and rearranging, we get |42| : 



O(e ) : f (0) = f e<? , 
Oie 1 ) : (d t0 + m)f (0) = -M (1) + S, 



0(e 2 ) : djW + idtv+m 



1 A 

2 



f(l) = _Af(2) 



(A55) 
(A56) 
(A57) 



where Ei = T(e ai I)T 1 , i G x,y. After substituting for f(°), Ei and S, the first-order 
moment equations, i.e. Eq. ( ]A56j) become 



d to p + d x (pu x ) + dy(puy) = 0, 



d to (pu x ) + d x ( -p + pu x ) + d y (pu x u yJ 



d tQ (pUy) + 8 X (pU X Uy) +0y [-P + PU 



(A58) 
(A59) 

(A60) 



^0 ( ^ + P( U l + U 2 y)^j + d X (J^PUX + PUxUl^ + Oy ^PUy + P^Uy 



-u 3 ff } + 2F x u x + 2F yU y, (A61) 



9t (p(u 2 x - u 2 y )) + d x {^pu x - pu x u 2 y ^j + d y y-T^pUy + pu 2 x u y ^j 



-ujf ] + 2F x u x - 2F y u 



(A62) 
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d tQ (pU X Uy) + 8 X -pUy + pU x Uy + Oyl - pU x + pU X U 



+ F x Uy + FyU X , 



(A63) 



9t {^P u y + P u l u y^j + Q» (p^%) + 3/ Qp + ^p( m * + wj) + P" 

= -w 6 /6 + F v u l + 2F x u x u y 



x u y 



(A64) 



<% Qpw* + + d x + ip(w^ + u 2 y ) + pulul^j 



d y (pu x u y 



'Mitt + F xU 2 y + 2FyUyU X , 



(A65) 



dt f~P + ~p(«* + u 2 y ) + pu x Uyj 



1 2 \ , o / 1 , .2 . 



+ I oP M ^ + PU X U y ) +dy[ -PUy + pU x Uy 



-ugff 1 + 2F x u x u 2 y + 2FyU y u x u 2 x . (A66) 



Similarly, the second-order moment equations can be derived from Eq. ( 1A57I) . which can 
be written as 

d to p = 0, (A67) 



dti (pux) + d x 
9tt (pu y ) + d x 
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7iX) 
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0, (A68) 
0, (A69) 



9 tl [ gP + + U y ) \ + 9 ifl 

+ d v 



+ 3* 
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d tl (pu x u y ) + <9 4o 



+ d v 



1 - ^5 ] if 



+ 9, 



1 

-c 
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1 " > ) F 



5 > 



(A72) 
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d h [ ^P u y + P u x u y j + 9 t 

+ d v 



1 - 2^6 ) /< 



TV) 



+ d x 



1 - 2^ ) ft 



5 



1 - 2 Ws ) fa 



70) 
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Qp + \p( u I + mJ) + pu^uj^ + d to 



1 - ^8 ] ^ 



1 " > ) ^ 



+ a, 



i 



1 - 2^6 ) /. 



H 1 ) 



(A75) 



Combining Eqs. (1A"58|) . (1A"59|) and ( 1A"60|) . with e times Eqs. ()A"67jl . (1A68]) and f lA"69|) . 

respectively, and using <9 t = <9 to + e<9 tl , we get the dynamical equations for the conserved or 
hydrodynamic moments after setting the parameter e to unity. That is, 



d t p + d x (pu x ) + d y (pu y ) = 



(A76) 
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+ F X , 



(A77) 



dt(p%) + d x (pu x u y ) + d y (pu 2 y ) = -d x (\p \ ~ 9 * f 1 ~ T^ 5 ) /{ 



'5 



<5„ 



1 



~ ( 1 - ~UJ 3 J ^ - =- ( 1 - > 4 J / 



+ F v . (A78) 



2 V 2 7" 2 V 2 
In the above three equations, Eqs. f)A76l) - flA78|) . we need the non-equilibrium raw moments 
fl 1 ^ and fl^ or tt X x + Kyy\ ^xx 1 — ^yy 1 and TTxy\ respectively. They can be obtained 
from Eqs. (1A62|) . ( IA63j) and (1A64j) . respectively. Thus, 



70) 

J3 



1 

U 3 



-dt Q ( 3P + P^X + Uy] 



Or ( ~PU X + pU X U 2 y \ - 8y \ \pUy + pU 2 x Uy 
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m 

5 



1 

-\-F x Uy ~\- FyU^ 



i 



-pW y + pW^ 



<9„ 



1 



(A81) 



The above three non-equilibrium moments can be simplified. In particular, by using the first- 
order hydrodynamic moment equations, Eqs. (1A58j) - (1A60l) and neglecting terms of 0(u 3 ) or 
higher, we have d to {pu 2 x ) « 2F x u x , d to (pu 2 y ) « 2F y u y and d to (pu x u y ) « + F^. Substi- 
tuting for these terms in Eqs. ( 1A79j) -( fA81j) . and representing the components of momentum 
for brevity as 

Jx PU X , Jy PUy, 

we get 





2 ^ -> 

3w 3 


(A82) 
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3w 5 ^ x3y + ' 
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Now, let 
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1 ( 1 



3 V w 5 2 



(A85) 



1 -- 1 

3 2/ 3 \o>4 2, 

and substituting the simplified expressions for the non-conserved moments, Eqs. ( 1A82I) - 
( lA84p . and by using the relations for relaxation parameters given in Eq. flA85]) in the con- 
served moment equations, Eqs. (lA76l) - flA78j) . we get 



(A86) 
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d t j y + d x 



Jxjy 
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+d v 



M^yjy-^ -7)+^ ■ 



V 



(A88) 



where p — ^p is the pressure field. Evidently, the relaxation parameters Co>4 and deter- 
mine the shear kinematic viscosity of the fluid, while U3 controls its bulk viscous behavior. 
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Moreover, w 4 = w 5 to maintain isotropy of the viscous stress tensor ($4 = t? 5 ). Thus, the pro- 
posed semi-implicit procedure for incorporating forcing term based on a specialized central 
moment lattice kinetic formulation is consistent with the weakly compressible Navier-Stokes 
equations without resulting in any spurious effects. 

It may be noted that in this work, we have employed a multiscale, or more specifically a 
two time scale, expansion j^j] to derive the macroscopic equations. An alternative approach 
is to consider a single time scale with an appropriate scaling relationship between space 
step and time step to recover specific type of fluid flow behavior. This broadly leads to 
two different types of consistency analysis techniques: (a) asymptotic analysis approach 46 1 
based on a diffusive or parabolic scaling [9] and (b) equivalent equation approach used in 



47|, 



48| based 



conjunction with certain smoothness assumption and Taylor series expansion 
on a convective or hyperbolic scaling 49]. A recursive application of the LBE and an 



associated Taylor series expansion without an explicit asymptotic relationship between the 
lattice parameters can also be used to analyze the structure of the truncation errors of the 
emergent macroscopic equations jsoj] . Another more recently developed approach is that 
based on a truncated Grad moment expansion using appropriate scaling with a recursive 
substitution procedure [36J , which has some features in common with an order of magnitude 
analysis for kinetic methods 51| . It is expected that such analysis tools can alternatively 
be applied to study the new computational approach described in this work. 



Appendix B: Generalization of Equilibrium and Sources with a Multiple Relaxation 
Time Cascaded Lattice Kinetic Formulation 

Let us first consider relaxation process of second-order non-conserved moments in the 
rest frame of reference: 

rp = ^9p = ^ a f ~(^ P \ = 3,4,5. (Bl) 

Here, summation of repeated indices with the subscript /3 on the RHS is not assumed and 
the superscript "c" for gp represents its evaluation for cascaded collision process, with g*p 
given in Eq. (IA43I) but restrict here to second-order moments. For convenience, we now 
define the non-equilibrium (raw) moment of order (m + n) as 
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or equivalently / 



fp — fp, where (3 = m + n. Thus, 



9z 



UJ3 
"12 



K xx ' K yy 



• 12 /3 , 



3S 



CU5 
' 4 



d(neq)' C^neq)' 



~J4 



w 5 9 (ne9) 

4 /5 , 



(B3) 
(B4) 
(B5) 



The next step is to relax the third and higher order non-conserved moments in the 
moving frame of reference, with each central moment relaxing with distinct relaxation time, 
in general. That is, 



E 





<(e, 



ax ^x 



,) m (e ay -u y ) n \Kp)gl = up 



K — K x myn + U 



x" l y" 



m + n>3. (B6) 



Clearly, this is equivalent to considering the last three rows of the J 7 matrix given in 
Eq. (1A41I) to determine g^, for = 6,7,8 [3^]. Expanding the terms within the brack- 
ets of the RHS Eq. (IB6I) in terms of raw moments, we get 



^xxy 


l^xxy 


@xxy 




~d(neq)' C:(ne?)' 
K xxy ~~ ^ U x K X y 


^_eq 

xyy 


^xyy 


— a xyy = 




"^(neq)' d(neq)' 
K xyy ~ ^ u y K xy 


^eq 

xxyy 


ft'xxyy 


~ @xxyy 




~d{neq)' d(neq)' 
^xxyy ^ U x K xyy 








+4:U x UyK xy 



u y K xx 

dneq) 1 ) 

- U X^yy 

_ 9 dneq)') 2 — ( ne lY) , 2— ( ne 9)') 
<ZU y K xxy + U x K yy U y li xx 



(B7) 
(B8) 

(B9) 



Now, in a manner analogous to the relaxation of second-order (raw) moments to their 
equilibrium states, we assume relaxation of third and higher order (raw) moments to their 
corresponding "equilibrium" states as well, which are as yet unknown, but will be determined 
in the following consideration. That is, we consider the ansatz 



,-jeq,G _ -j , . 

% = up a ,„ M , P = 6,7,8. 



{Kp\Kp) 



(BIO) 



Here, the superscript "G" represents the "generalized" expression, i.e. f a ' is the generalized 
equilibrium in the presence of forcing terms (due to the presence of the 'overbar' symbol), 
which for a = 6, 7, 8 will be determined in the following. Again, summation of repeated 
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indices with the subscript j3 on the RHS is not assumed. Evaluating Eq. (1B10I) yields 

^eq,G' 



4 

97 = — 



xxy xxy 



~eq,G' d 
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■>:yy 
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UJ8 

4 L 



K xyy 

~eq,G' 
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^eq 
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(Bll) 
(B12) 
(B13) 



Now substituting Eqs. (1B2I) . (1B7I) -( 1B9I) and ( 1B10I) in Eq. (1B6I) and simplifying and rearranging 
the resulting expressions yield the desired expressions for the generalized equilibrium in the 
presence of forcing terms 



^eq,G' 
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~eq,G' 
K xyy 

~eq,G' 
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(B14) 
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(B16) 



where the coefficients ^ a in Eqs. (IB14[) - flB16[) are functions of the various ratios of the 
relaxation times of the above MRT cascaded formalism and velocity field arising relaxing 

eq,G' 



the moments in the moving frame of reference. The coefficients for « are 



4 

¥>6 



and tor k„,„, are 
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^ (1 " *S) «y, 
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and, finally, for K^xyy are 
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3 + + V 
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^8 



y 



u x . 



2 (i - el) 

Here, in Eqs. (lB17j) - flB19l) . the parameter 02 refers to the ratio of relaxation times u a and 
u)p. That is 

(B20) 



Up 
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_ieq,G' — G 

Now, in the notations of the previous section, we can rewrite K xmyn in terms of f ^ , or 



more explicitly, in terms of the regular generalized equilibrium and source moments, i.e. / G 



and Sp, respectively, using f ^ = fp — \S G ' . Thus, compactly, the generalized equilibrium 
and source moments are 



N v N v 



Q=3 

N v 



a=3 



a=3 



where N„ 



S° = S P -J2^S«, = 6,7 
5, = 6,7 



(B21) 
(B22) 



7, = 8 



. It should, however, be noted that f e g q,G = fl q and S9 = Sp for 



(3 < 5. This analysis further extends that of Asinari [36[], who showed generalized equilib- 
rium for a particular form of Cascaded-LBM without forcing terms. Thus, the generalized 
equilibrium arising from the cascaded nature of the collision step for the third and higher 
order (raw) moments is a function of conserved moments, non-equilibrium part of the lower 
order moments and the various ratios of the relaxation times in the MRT formulation. Sim- 
ilarly, the generalized sources for the third and higher order moments is a function of the 
products of force fields and fluid velocity, as well as the ratio of relaxation times. In view of 
the above, the cascaded formulation can also be reinterpreted by defining the generalization 
of the equilibrium and source in terms of the following local coefficient matrix C = C(lt, t): 
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That is, if the information cascades from lower to higher moments during a time interval 
(t, t + 1), the raw equilibrium and source moments in the lattice frame of reference generalize 
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to 



= (l-W^+Cf^D, (B24) 



= (X - C) Srt A (B25) 



where t* represents some intermediate time in (t, t + 1). Clearly, the generalization of both 
equilibrium and sources degenerate to corresponding regular forms only when the relaxation 
times of all the moments are the same. That is, when the approach is reduced to the SRT 
formulation, f^ g ' G = f^ q and = Sp for all possible values of /3, since C = 0, i.e. a null 
matrix in that case. In the previous section, a consistency analysis for a special case of the 
central moment method was presented. The same notation and procedure can be adopted for 
the general case involving cascaded relaxation (represented as a relaxation of non-conserved 
raw moments to their generalized equilibrium) with generalized sources presented here, when 
fp q becomes f^ 9 ' G and (l — Sp becomes Sp for f3 = 6,7,8. Inspection of the details 
of the Chapman- Enskog moment expansion analysis presented in the earlier section shows 
that the consistency of the Cascaded-LBM to the NSE remains unaffected by the presence of 
generalized equilibrium and sources. In particular, the generalized forms contain coefficients 
which are functions of local fluid velocity and the ratio of various relaxation times, and 
terms that are non-equilibrium part of the lower order moments, which are negligibly small 
in nature for slow or weakly compressible flows, as they involve products of various powers 
of hydrodynamic fields. Since for consistency purpose, we need to retain only 0(Ma 2 ), the 
presence of the generalized terms do not affect the end result of the derivation presented in 
the previous section. 

An interesting viewpoint to note is that the use of relaxation to generalized equilibrium 
(including the effect of sources), i.e. Eq. (IBlOj) may be considered as an alternative com- 



putational framework to actually execute the cascaded MRT collision step. It reduces to a 
corresponding TRT collision step, when o; even = Co> 4 = Uq = u 8 and u odd = u 3 = u 5 = u 7 . 
Also, a different perspective of the generalized equilibrium, Eq. ( 1B21I) can be arrived at in 



light of the consistency analysis performed in the previous section. For example, for the 
third-order moments, /3 = 6 and 7, Eq. (1B21|) needs the non-equilibrium moments f^ neq \ 
ft eq) and fl neq \ which can be approximated by Eqs. flA82j) . flA83j) and (IA84j) . respectively, 



which actually provide expressions for the components of the strain rate tensor in the cas- 
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caded formulation. Thus, we get 



f!^ « ft ~ ~ Pi Uy^ ■ 7 ~ - Pj Uy (d xJy - d y3x ) 

u x {d x j y + d y j x ) , (B26) 



3 V^5 w 6 

f^ G « ff q - - ( ) Ux V ■ 3 - £ ( 1 u x (d x j y - d y j x 





r- 


1 




V^4 


W 7 



-of- - ) M 2/ (Sail/ + fyjas) • ( B27 ) 

In other words, the generalized equilibrium is a function of density and velocity fields and 
their gradients, the coefficients of the latter terms are given as difference of relaxation times 
of moments of different order. 



Appendix C: Introducing Time-implicitness in the Cascaded Collision Operator 

Here, let us investigate the possibility of developing an executable LBE formulation where 
implicitness in time is introduced in the cascaded collision kernel, which could be useful in 
certain applications. In particular, we extend Eq. (|35|) such that the cascaded collision 
operator Fl^nt ^ is now treated to be semi-implicit in time: 

f a (~£ + ~t a , t + 1) = f a (at, t) + - [{K ■ g) a (x>,t) + (£ • B)aO*+-£„,*+l)] 

In order to avoid an iterative procedure for the use of Eq. ( 1C1I) . we now define the following 
transformation with the introduction of a new variable h a : 

h a = f a -^(K-E) a -^S a . (C2) 

Now, substituting Eq. flC2]) in Eq. ([CIj> . we get 

h a (~$ + lt a , t + 1) - hjj£, *) = (£• g) a &,t) + S a & tt ) (C3) 

As a result, Eq. (1C3|) now becomes effectively explicit. In the new variable, the hydrodynamic 
fields can be obtained as p = Yla=o h a and put = Yla=o h a e a i + \Fi. The post-collision values, 
i.e. h a can be obtained by replacing f a with h a in Eqs. (jl65p -( JTT3"|) . Now, to obtain the 
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collision kernel {K ■ g) a in Eq. (1C3I) in terms of we define the following raw moment of 
order (m + n): 

T] x m y n = ^ ] h a e ax e a y = (& a x e a.y\ha) j (^4) 

Q 

where can be represented and computed in a manner similar to that given in Eqs. fll2ip ~ 
ffT26D . From Eqs. (|C2]) and (JClJ), we obtain 

Vx m y n = K x m y n ~ 7} ^ ] (-^/3 1 e ax e ay) 9fi ~ q ". 



= %™ y n (KpKxCy) 9fi (C5) 

where (^le^e^) ^3 can be obtained by exploiting the orthogonal properties of /C, i.e. 
from Eqs. (155]) -(155]). 

Now substituting Eq. (IC5p in the collision kernel written in compact notation as given in 
Appendix \A\ i.e. in Eqs. ( 1A31j) - (1A36j) . and simplifying we get 



1 CU3 J 2 / 2 2 \ /— ; — ' \ N ! 

#3 = 77: 7^ I 1 V i oP + P{U X + UJ - [Vxx + Vw) ~ n(°xx + O \ > (C6) 



12(1 + 1^3) \?> KV * 2/7 Wxa; lyy ' 2 

^4 f 2 2N /S' S' N l . 



^ = 4 (1 + 1^) ~ M ^ ~ ~ ^ ~ 2^ xx ~ ° yy) j ' (C7) 

% = i (ifr^y { - ty - \* xy ) , (C8) 

#6 = ^ ^ _|_ l y ^ ^pM^My + 77 x . x . y — 2u x 1] xy — U y T] xx — 2^ 



~0%( 3 ?3 + &) - 2^5, (C9) 



^ 1 0J 7 ( 2 - - - 1- 

-\^x{Zgz - g A ) - 2u y g 5 , (CIO) 



1 w 8 fl . q 2 2 

1 + 3pu x u y - 



Vxxyy ^' U x r lxyy < ^ U yVxxy U xVyy U yV. 



4(l + |o; 8 ) 1 9' 
+4u x u y rj xy 

-4u x u y g 5 - 2u y g 6 - 2u x g 7 . (Cll) 



It may be noted that a Chapman-Enskog analysis, as given in Appendix \A\ when per- 
formed with the above collision operator, yields the following relations between relaxation 
parameters and transport coefficients (see Eq. (1A85I) ): $3 = $ 4 = $5 = 
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for the hydrodynamical equations given in Eq. ( 1A87j) and ( 1A88j) . Thus, the above con- 
siderations show that it is possible to introduce time-implicitness in the cascaded collision 
kernel, and when a transformation is introduced to make the computational procedure effec- 
tively explicit, it leaves the form of "gp unchanged with a simple re-scaling of the relaxation 
parameters. 
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